NASA/TM-20 14-2 18175 



First-Order Hyperbolic System 
Method for Time-Dependent 
Advection-Diffusion Problems 


Alireza Mazaheri 

Langley Research Center, Hampton, Virginia 
Hiroaki Nishikawa 

National Institute of Aerospace, Hampton, Virginia 


March 2014 



NASA STI Program ... in Profile 


Since its founding, NASA has been dedicated to the 
advancement of aeronautics and space science. The 
NASA scientific and technical information (STI) 
program plays a key part in helping NASA maintain 
this important role. 

The NASA STI program operates under the 
auspices of the Agency Chief Information Officer. 

It collects, organizes, provides for archiving, and 
disseminates NASA’s STI. The NASA STI 
program provides access to the NASA Aeronautics 
and Space Database and its public interface, the 
NASA Technical Report Server, thus providing one 
of the largest collections of aeronautical and space 
science STI in the world. Results are published in 
both non-NASA channels and by NASA in the 
NASA STI Report Series, which includes the 
following report types: 

• TECHNICAL PUBLICATION. Reports of 
completed research or a major significant phase 
of research that present the results of NASA 
Programs and include extensive data or 
theoretical analysis. Includes compilations of 
significant scientific and technical data and 
information deemed to be of continuing 
reference value. NASA counterpart of peer- 
reviewed formal professional papers, but 
having less stringent limitations on manuscript 
length and extent of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 
and technical findings that are preliminary or of 
specialized interest, e.g., quick release reports, 
working papers, and bibliographies that contain 
minimal annotation. Does not contain extensive 
analysis. 

• CONTRACTOR REPORT. Scientific and 
technical findings by NASA-sponsored 
contractors and grantees. 


• CONFERENCE PUBLICATION. 

Collected papers from scientific and 
technical conferences, symposia, seminars, 
or other meetings sponsored or co- 
sponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 
technical, or historical information from 
NASA programs, projects, and missions, 
often concerned with subjects having 
substantial public interest. 

• TECHNICAL TRANSLATION. 
English-language translations of foreign 
scientific and technical material pertinent to 
NASA’s mission. 

Specialized services also include organizing 
and publishing research results, distributing 
specialized research announcements and feeds, 
providing information desk and personal search 
support, and enabling data exchange services. 

For more information about the NASA STI 
program, see the following: 

• Access the NASA STI program home page 
at http://www.sti.nasa.sov 

• E-mail your question to help@sti.nasa.gov 

• Fax your question to the NASA STI 
Information Desk at 443-757-5803 

• Phone the NASA STI Information Desk at 
443-757-5802 

• Write to: 

STI Information Desk 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320 


NASA/TM-20 14-2 18175 



First-Order Hyperbolic System 
Method for Time-Dependent 
Advection-Diffusion Problems 


Alireza Mazaheri 

Langley Research Center, Hampton, Virginia 
Hiroaki Nishikawa 

National Institute of Aerospace, Hampton, Virginia 


National Aeronautics and 
Space Administration 

Langley Research Center 
Hampton, Virbinia 23681-2199 


March 2014 



The use of trademarks or names of manufacturers in this report is for accurate reporting and does not 
constitute an official endorsement, either expressed or implied, of such products or manufacturers by the 
National Aeronautics and Space Administration. 


Available from: 

NASA Center for AeroSpace Information 
7115 Standard Drive 
Hanover, MD 21076-1320 
443-757-5802 


Abstract 


A time-dependent extension of the first-order hyperbolic system method [J. 
Comput. Phys., 227 (2007) 315-352] for advection-diffusion problems is in- 
troduced. Diffusive/viscous terms are written and discretized as a hyperbolic 
system, which recovers the original equation in the steady state. The resulting 
scheme offers advantages over traditional schemes: a dramatic simplification in 
the discretization, high-order accuracy in the solution gradients, and orders-of- 
magnitude convergence acceleration. The hyperbolic advection-diffusion sys- 
tem is discretized by the second-order upwind residual-distribution scheme in 
a unified manner, and the system of implicit-residual-equations is solved by 
Newton’s method over every physical time step. The numerical results are 
presented for linear and nonlinear advection-diffusion problems, demonstrat- 
ing solutions and gradients produced to the same order of accuracy, with rapid 
convergence over each physical time step, typically less than five Newton iter- 
ations. 
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1 Introduction 


In this paper, we introduce an unsteady hyperbolic advection-diffusion scheme. 
The first-order hyperbolic system method was initially proposed in Ref. [1] 
as a radical approach to steady diffusion problems: discretize an equivalent 
first-order hyperbolic system to solve the diffusion equation. It has a number 
of attractive features such as the discretization of high-order derivatives by 
methods for hyperbolic systems, equal order of accuracy for the solution and 
the gradients ( viscous /heat fluxes), and orders of magnitude acceleration in 
steady convergence. These advantages have been demonstrated for the diffu- 
sion equation [1] and the advection-diffusion equation [2] by the second-order 
residual-distribution method, and for the compressible Navier-Stokes equa- 
tions [3] by the second-order finite-volume method, and for the advection- 
diffusion equation by first, second and third order finite- volume methods [4,5]. 
These lines of development, however, have been restricted exclusively to steady 
state problems. Towards enabling accurate practical time-dependent computa- 
tions by the first-order hyperbolic system method, this paper presents the first 
study on the unsteady extension of the first-order hyperbolic system method. 
We demonstrate that time-accurate schemes constructed via the hyperbolic 
method produce high-order accurate solutions and derivatives at every time 
step and allow the construction of a highly efficient inner solver for implicit 
time stepping. 

Time-dependent computations are possible by implicit time-integration 
methods such as the second-order backward-difference scheme, which is widely 
used for practical computations. To perform the implicit time integration, it 
is required to solve a system of global time-dependent residual equations over 
each time step. We employ a steady solver developed by the first-order hy- 
perbolic system method to solve the system efficiently and accurately. The 
spatial discretization is performed by the residual-distribution method as in 
Refs. [1,2]. Taking advantage of the compactness of the residual-distribution 
schemes, we develop a fully- implicit steady solver with exact linearization, i.e., 
Newton’s method, for a second-order upwind scheme. The strategy taken here 
may be thought of as a dual-time formulation [6]. However, we introduce 
a pseudo-time just for convenience of discretization: the first-order system 
is hyperbolic in the pseudo-time and thus will be discretized by the upwind 
residual-distribution scheme. Once the spatial discretization is complete, the 
pseudo-time step is taken to be infinity to define the system of time-dependent 
residual equations, which is then solved by Newton’s method. In this paper, 
we consider one-dimensional linear and nonlinear advection-diffusion equa- 
tions. These examples are sufficient to illustrate the general methodology 
for enabling time-dependent computations by the hyperbolic method; this is 
the main goal of this paper. The introduction of this unsteady hyperbolic 
advection-diffusion scheme along with its great advantages is served as a foun- 
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dation for the extension of the scheme to higher-dimensions and more complex 
systems. 

It is emphasized also that there are real practical applications that one- 
dimensional analyses are sufficient and routinely being used in industries (e.g., 
NASA). For instance, in-depth material thermal response calculations of ther- 
mal protection systems of atmospheric entry vehicles, where advection-diffusion 
equation is solved through the material, is a prime example of such applica- 
tions [7-9]. Another example is the experimental aeroheating data reduction 
that are obtained with global phosphor thermography techniques [10,11]. Ef- 
ficient and accurate one-dimensional schemes as presented in this paper will 
significantly improve these analyses and could potentially make such calcula- 
tions part of the routine aerothermodynamic environment predictions. The 
present paper serves, therefore, also as an important foundation for the exten- 
sion to more practical one- dimensional problems. 

There exist other methods for enabling time-accurate computations by 
the residual-distribution method: the Crank-Nicolson method [12], the space- 
time method [13-15], the backward-difference method [16,17], and the explicit 
Runge-Kutta-type method [18]. These methods are equally applicable to the 
construction of time-accurate schemes in the hyperbolic method. In this pa- 
per, we employ the backward-difference method for simplicity but without los- 
ing practical applicability. The main difference between our time-integration 
method and those in Refs. [16, 17] lies in the method for solving the sys- 
tem of residual equations arising from the implicit time integration scheme: 
we employ Newton’s method whereas Ref. [16] employs the dual-time step- 
ping method with a pseudo-time marching iteration and Ref. [17] employs the 
multigrid method with a point-implicit iteration. 

The paper is organized as follows. In the next section, the hyperbolic 
advection-diffusion system is described. In Section 3, a compact second-order 
residual-distribution scheme, a steady solver, and the boundary conditions are 
discussed. In Section 4, the construction of time-accurate schemes is given. 
In Section 5, the time-accurate scheme is extended to a nonlinear advection- 
diffusion equation. In Section 6, numerical results are presented. Finally, 
Section 7 concludes the study with remarks. 


2 Hyperbolic Advection-Diffusion System 

Consider the one-dimensional (1-D) advection-diffusion equation, 

d t u + a d x u = v d xx u , (1) 

where a and v are both taken to be positive constant. We will follow the 
procedure described in Ref. [2] and rewrite the above equation with a first- 


4 



order advection-diffusion system: 


( 2 ) 

( 3 ) 


d t u = —a d x u + v d x p, 
d t p = (d x u - p)/T r , 

where the relaxation time, T r > 0, is arbitrary. Towards the steady state, 
the variable p approaches the solution gradient and hence the above equation 
becomes identical to the original advection-diffusion equation in the steady 
state. In the vector form, the first-order advection-diffusion system can be 
written as 

d t U + Ad x V = S, (4) 

where 



a —v 

— 1 /T r 0 J ’ 


S 


0 

- P/T r 


( 5 ) 


The above system is hyperbolic in time because it has the following real eigen- 
values for any positive T r : 



with linearly independent eigenvectors [1]. 

Note that the above hyperbolic formulation is equivalent to the original 
Eq. (1) only in the steady state. The time derivatives in the hyperbolic sys- 
tem, at least d t p, should therefore be taken as pseudo-time derivatives. They 
serve mainly as a guide for discretization by making the whole system hy- 
perbolic in time for which various discretization techniques are available, e.g., 
upwinding. The benefits are, however, much more than just the convenience in 
discretization. First, the hyperbolic discretization typically results in a strong 
coupling between the two variables and achieves the equal order of accuracy 
for u and p on arbitrary grids. Second, the resulting explicit numerical scheme 
is stable with an 0(h) time step through the diffusion limit or 0(1/ h) condi- 
tion number in the residual Jacobian for implicit solvers. It has been shown 
to yield 0(1 /h) acceleration in convergence over traditional methods for the 
diffusion equation [1,4], for the advection-diffusion equation [2,5], and for the 
compressible Navier-Stokes equations [3]. Unsteady computations are possible 
by implicit time stepping, which is the main subject of this paper. 

In the next section, we first fully describe the advection-diffusion equation 
discretization process, the implicit steady state formulation, and the boundary 
condition implementation in the steady state limit. We then extend the scheme 
to time-accurate and nonlinear advection-diffusion equation in the sections 
that follow. 
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Figure 1: Schematic of grid spacing for a 1-D grid. 


3 Discretization, Steady Solver, and Bound- 
ary Condition 

3.1 Discretization 


To discretize the first-order hyperbolic advection-diffusion system, we employ 
Residual-Distribution (RD) method. The method consists of two steps; 1) 
residual evaluation over the cells (or elements), and 2) distribution of the 
residuals to the nodes. 

Consider a one-dimensional domain discretized with N nodes that are dis- 
tributed arbitrarily over the domain of interest with the solution, u, and the 
solution gradient, p, data stored at the nodes denoted by Xi, i — 1,2, 3, ..., N 
(Fig.l). We define the cell-residual <F C by integrating the spatial part of the 
Eq. (4) over the cell, C , defined by the nodes i and i + 1: 


PXi+l 


—a(u i+ 1 - Ui) + z/(p i+1 - Pi) 


<p c = 


-A<9 X U + S)dx = 


T 
± r 


Pi + 1 +Pi/ \ 

^2+1 ^ i 0 y^i- 1-1 %i) 


P) 

Note that the source term integration has been evaluated by the trapezoidal 
rule to ensure second-order accuracy. It should be noted also that the flux bal- 
ance terms e.g., u x , has been integrated exactly, which is not possible in higher 
dimensions and a high-order quadrature would be required to achieve the equal 
order of accuracy for the solution and the gradients in the residual-distribution 
method as described in Ref. [2]. We again remark that the relaxation time, 
T r , is arbitrary because in the steady state limit p = u x for any T r . We com- 
plete the evaluation of the cell residual with the definition of T r . Using the 
dimensional analysis [2], this is defined as 


T 
± r 


max(Ai, A 2 ) 


a 

2 




( 8 ) 


where L r is a length scale that may be optimized using Fourier analysis to 
enhance the convergence [19]. The effect of the optimized length scale remains 
to be investigated. Here we use the value recommended in Ref. [1,4,5]: 


/. . — 


1 

27t' 


( 9 ) 
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We now solve for T r from Eq. (8): 

T = 

1 r 


( 10 ) 


a + is/L r 

Substituting back the relaxation time into Eq. (6), we express the eigenvalues 


as 


\i = —v/L r , A 2 — a + u/L r . (11) 

At this point the residual evaluation is completed. We now define the 
distribution of the residuals. We achieve this by first diagonalizing the matrix 
A with the following right and left eigenvectors: 


R = 


— Ai T r — A 2 T r 

^ aL r is 

1 —L r 

1 1 

L r (aL r H - 2is) 

— 1 (uL r )/(aL r + u) 


.( 12 ) 


We then expand the matrix A to two separate parts that are distinguished by 
its corresponding wave speeds (eigenvalues): 


where 


A = RAL = Ai B~ + \ 2 B + , 


A 


Ai 0 
0 A 2 ’ 


(13) 

(14) 


BT 


1 

aL r -\- 2is 


is isL r 

(aL r “h is) j L r aL r -f~ is 


(15) 



aL r — | — 2 zv 


aL r + v —vL r 
— (aL r + u)/L r v 


(16) 


The matrices B + and B~ project the solution, respectively, onto the left- 
and right-running waves. Hence, this is upwinding. We also note that the 
distribution matrices have the following property 


B~ + B + = /, 


( 17 ) 


which is required for conservation. We can now distribute the cell residual 
3> c to the nodes using the projection matrices B + and B~ as described in 
Fig. 2. The distribution is done this way because the left running wave, A 2 , 
for example, is contributing to the B + and thus the cell residual on the left cell 
is weighted with the B + . After the distribution step, we obtain the following 
semi-discrete scheme in pseudo-time: 

= 1(H+$ l + H-$ r ), (18) 

dr hi 

where <h L and Q R denote the cell-residuals over the left and right cells of the 
node i, respectively, and h % is the dual volume (see Fig. 1) defined by 


hi 


hL + hn 


h'L 1) hji 2)j+l 


(19) 
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Figure 2: Residual distribution to the nodes. 


3.2 Steady Solver 

Equation (18) can be solved with either an explicit or an implicit formula- 
tion. We choose to focus on the implicit formulation for a better transition 
to the next sections, where we extend this scheme to general nonlinear time- 
dependent advection-diffusion problems. Therefore, we drop the pseudo-time 
derivative from (18) and solve the resulting system of steady residual equa- 
tions, 


0 = Res, 


( 20 ) 


by a fully implicit solver. We remark that in some complex problems the 
pseudo time derivative should be kept for stability purposes. In our one- 
dimensional problems, however, we could choose an infinite pseudo time step 
without encountering any stability issue. 

The implicit formulation is defined by 

U fc+1 = U fc + AU fc , (21) 

where U = (u\,pi,U2,P2, ■ ■ ■ ,Un,Pn) and k is the iteration counter. The 
correction AU fc = U fc+1 - U k is determined as the solution to the linear 
system: 


= -Res‘, (22) 

where Res^ is the steady residual vector evaluated by U fc . The Jacobian 
matrix is sparse, having three 2x2 blocks in each row for all interior nodes 
and two blocks for boundary nodes. The z-th interior pair of row of the linear 
system is given by 


ARt, + J,AU k + .J l+i AU k +l = ~Ub + O l + B~O r )\ (23) 

lii 


where A Uf_ x = Ap£_i), A Uf = (Av$,A$), AUf + i = (AuJ +1 , Apf +1 ), ; 


Ji—l 

Jr 

Ji+1 



d$ L 

dU~i 


1 

hi 



d$ L 

m 


+ b- 


d$ R \ 

~m) ’ 


1 

hi dU i+ i 


and the derivatives of the cell-residuals are given by 


d$ R 

~dUi 




(24) 

(25) 

(26) 

(27) 


T —a v 

d$ L _ 

Wl ~ JL h L 

T r 2 T r 


(28) 


The linear system may be solved efficiently by Thomas’ algorithm, which 
is an 0(N ) method. In this paper, however, we employ the Gauss-Seidel (GS) 
relaxation, which is also an O(N) method for the discretization arising from 
hyperbolic advection-diffnsion system and extends straightforwardly to more 
complex systems, irregular grids, and higher dimensions, whereas Thomas’ 
algorithm is not applicable in higher dimensions. The GS relaxation scheme 
may be applied equation by equation (decoupled relaxation) or node by node 
(collective relaxation) updating the solution set at a node simultaneonously 
[20]. In this work, the collective GS relaxation is employed for robustness; the 
decoupled GS relaxation was found to give better convergence in most cases 
but also found to be unstable in some cases whereas the collective relaxation 
encountered no such problems. In actual computations, we do not fully solve 
but relax the linear system to reduce the linear residual by two or three orders 
of magnitude. 

Note that the Jacobian is exact and thus the implicit algorithm is New- 
ton’s method. It is one of the advantages of the residual-distribution method 
that second-order accuracy is achieved within a compact stencil, thus allowing 
the exact linearization with a sparse Jacobian matrix. The same is true for 
extension to higher dimensions. For linear problems, the method is essentially 
a direct method. That is, if we solve the linear system exactly, then we get 
the solution at the first Newton iteration. 
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Figure 3: Schematic of the Dirichlet and the Neumann boundary conditions 
implementation in the hyperbolic advection-diffusion scheme. 


3.3 Boundary Condition 

In this section we describe the boundary condition implementation for our first- 
order hyperbolic system scheme. Because in our hyperbolic advection-diffusion 
system the second equation in the hyperbolic system solves for the solution 
gradient, p, the Neumann boundary condition is implemented exactly in the 
same way as the Dirichlet boundary condition. Therefore, we only show the 
formulation for the u variable and the same can be repeated for the p variable. 
Note that the discrete problem has a unique solution for linear problems with 
two values given at the left and right boundaries (this is consistent with the 
differential equation allowing two boundary conditions): 2(iV— 1) cell-residuals 
for 2N — 2 unknowns, which is the case in both of the following boundary 
conditions. 


Dirichlet or Neumann Boundary Condition 


For Dirichlet or Neumann boundary condition type schematically shown in 
Fig. 3, the discretized equation for the boundary condition imposed on u takes 
the following form at the first node: 

JfAC/f + J* AD 2 fc = (29) 

hi 

where for the imposed boundary condition on the u variable, the Jacobian 
matrices are 


j r 

r 2 


i B ~ 


l B ~ 


1 0 

d$ R (2) d$ R (2) ’ 

du\ dpi 

0 0 

d$ R (2) d$ R {2) ’ 

du 2 dp 2 


(30) 


(31) 


where <F R (2) denotes the second component of <& R corresponding to the second 
equation (the p equation) of the first-order hyperbolic system. The same 
process is repeated for imposed boundary condition on the last node (for either 
u or p.) 
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Figure 4: Schematic of the periodic boundary condition implementation in the 
hyperbolic advection-diffusion scheme. 


Periodic Boundary Condition 

To implement a periodic boundary condition we set Utv = Ui in Eq. (23) and 
obtain the following discretized equation and applied it on both two boundaries 
(Fig. 4): 


JN-iAU*^ + J,AE/f + J 2 AU\ = J B+ * L + B (32) 

where L denotes the cell defined by the nodes N — 1 and N, and R denotes 
the cell defined by the nodes 1 and 2. 


4 Extension to Time-Dependent Problems 

We extend the steady state formulation by modifying the first-order hyperbolic 
system given by Eq. (2) and rewriting it as a dual-time step system: 

d T u = —au x + vd x p — — u + s, (33) 

drP = ( d x u-p)/T r , (34) 


where t is the physical time, At is the physical time step, r is the pseudo-time, 
and s and a both depend on the physical time-step discretization scheme. 
Adopting Erst- or second-order Backward-Differencing-Formulation (BDF1 or 
BDF2, respectively) for the physical time integration, the a and s(x) terms 
are: 

u n 

a — 1, s = — : BDF 1, 


_ 2At n ~ 1 + A t n _ AC" 1 + At n n At"" 1 w _! 

° - At n - 1 + At n ’ 3 ~ At n - 1 At n U At n (At n - 1 + At n ) U ' 3 

(35) 

where a nonuniform time step formulation is used in the BDF2 formulation. 
The BDF1 is only used at the very beginning of the simulation with an initial 
time interval of A t°(« A t n ). After the first time interval, BDF2 with a 
uniform time step, At, is used. This procedure will ensure a second-order 
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accurate result through all times. The use of BDF1 is necessary because 
explicit time stepping is not possible with the hyperbolic system method. 

We solve the above system over each physical time interval to obtain the 
steady state solution in the pseudo-time space. In the steady state in r, we 
have p = d x u and therefore recover the consistent implicit time discretization 
scheme for the original time-dependent advection-diffusion equation. 

In the dual-time step system, the physical time derivative term acts as a 
source term to the first equation. This system thus shares the same eigenvalues 
and eigenvectors as the steady hyperbolic advection-diffusion system. We will 
see this more clearly by writing the time-dependent hyperbolic system in the 
vector form: 


<9U <9U 

d t dx 


DU + S, 


(36) 


where 


A 


j a —v 

, D = 

—a/ At 

0 

C - 

s(x) 

l 

1 

-J J 

o 

0 

-1/T r 

j ° — 

0 


(37) 


Because of the similarity in wave definitions in both time-dependent and steady 
hyperbolic systems, the distribution matrices B~ and B + are the same for both 
of these systems. Once the discretization is obtained, we drop the pseudo-time 
derivatives and solve the resulting system of discrete equations by Newton’s 
method as described in Section 3.2. To develop a time-dependent scheme, we 
therefore only need to add the source term effect in the cell-residual. 

The cell-residual <J> C of the time-dependent hyperbolic advection-diffusion 
system is 


= 


rxi+l 


-AU. C + DU + S )dx 


a 


- a(u i+ i - m ) + v{Pi+i - Pi) - (xi+i - Xi) — (u i+ 1 + Ui)/ 2 


fc+i 


T 

± r 


Pi+l + Pi ( 

^i + 1 2 (*^i+l 


+ 


(x i+1 -Xi)(s i+1 + Si)/2 


0 


(38) 


where k and n are the Newton iteration counter and the physical time index, 
respectively. Note that the second term of the Eq. (38), which is computed at 
the two previous physical time steps, is constant during the Newton iteration, 
and thus will not contribute to the Jacobian. The residual Jacobians needed 
in the Newton solver are exactly in the same form as Eq. (26), but the deriva- 
tives of the cell-residuals now include the contribution from the physical time 
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derivative: 


a 


—v 


d$ R 

dUi 


d$ L 

dUi 


a — h r 


2 At 


1 

± r 

—a — hj j 


a 

2At 


1 

T 7 

1 r 


h R ’ 
2 T r . 

v 

~2T r . 


(39) 


(40) 


The hyperbolic advection-diffusion scheme is now time-accurate. At each 
physical time level n, we solve the pseudo steady problem by Newton’s method 
with the current solution at n as the initial solution. 


5 Extension to Nonlinear Advection-Diffusion 
System 

Consider the generalized time-dependent nonlinear advection-diffusion equa- 
tion: 

d t u + d x f = d x (vu x ), (41) 

where / is a nonlinear function of u and v = v{u) is a function of u. We write 
the above equation using the dual-time stepping formulation in the following 


first-order hyperbolic system as 


d T u = -d x f + d x p - -^u + s, 

(42) 

T r 

, ,d T p = d x u p/v(u). 
u(u) 

(43) 

This is a preconditioned conservative form introduced in Ref. [3] to extend the 
hyperbolic method to nonlinear equations. Note that the variable p will be 

the diffusive flux in the pseudo steady state, not the gradient, 
form, the preconditioned conservative system is written as 

In the vector 

P -1 U r = — F x + S, 

(44) 


where 

P 1 


1 0 
0 T r /v(u ) 


F 


f~p\ 

-u )' 


S 


-£u + s(x) \ 

-p/v{u) ) ’ 


(45) 


and P is the preconditioning matrix. The main role of the preconditioning 
matrix here is not to optimize the condition number of the differential system, 
but rather to simplify the discretization. In particular, in this form, there is 
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no need to differentiate the diffusion coefficient with respect to the solution to 
derive the flux Jacobian. Multiplying both sides by P from the left, the flux 
term becomes 

OF. 


PF X = P^U X = AU„ 


where 


A = 


du 

a 

—v/T, 


-1 
r 0 


(46) 


(47) 


Here, we have introduced a = df /du, which is not a constant but a func- 
tion of u. Note that the preconditioned Jacobian is formally equivalent to 
the Jacobian of the linear time dependent hyperbolic ad vection- diffusion sys- 
tem. Hence, the eigenvalues and eigenvectors of this system are identical to 
those of the linear equations. It greatly simplifies the construction of the dis- 
tribution matrices as described below. This is one of the advantages of the 
preconditioned conservative formulation. 

To define the cell- residual, we first integrate the right hand side of Eq. (44) 
to get the unpreconditioned residual : 


tfC = 


r%i+ 1 


[— F-j, + S )dx 


' Xi 


= — (F, :+1 - Fj) + S<+1 2 +Si (s i+1 - X,). (48) 

For the advective flux that is quadratic in the solution, it is well known that 
a conservative linearization exists: 


* C = -||j(U. + i - U.) + S ‘ +1 2 + S ‘ (* 1+1 - x,), (49) 


where IP) is the analytical Jacobian evaluated at the arithmetic average of the 


solution, (Ui + 1 + Uj )/ 2 and it satishes 


dF 


Fj+l — Fj — ^0r(U i+ i — Uj), 


(50) 


exactly. The preconditioned cell-residual, which is to be distributed, is then 
given by 




c 


(51) 


where the matrix P is evaluated at the arithmetic average of the solution so 
that we have 


4> c = -A(U i+1 - UO + S?:+1 0 + (x i+1 - Xi ). 


(52) 
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The distribution matrices can then be constructed in the same way as in the 
linear case based on A: 


Bn = 


c aL r T 2v 


ULr 


(aL r T is) I L r aL r T is 


(53) 



1 

aL r 2is 


aL r + is —isL r 
— (aL r + v)/ L r V 


(54) 


where 


a = 


+ Q'i+1 

2 


IS = 


Vi + Vi + 1 
2 


(55) 


The subscript C indicates that these matrices are defined over the cell C. 
Note that the distribution matrices are not globally constant in general for 
nonlinear equations. The construction of the distribution matrices based on 
the conservative linearization is important for proper upwinding. If different 
averages are used, the distribution may not be strictly upwind. If the target 
equation does not allow the conservative linearization, the construction pro- 
posed in Ref. [21] may be employed to design a proper upwind distribution 
matrix. Note also that we evaluate is in the source term by the average value 
over the cell to simplify the second component of the cell-residual as follows: 


$ 


c 


(X ~\ 

-(/*+ 1 - fi ) + \Pi + 1 - Pi) ~ ( x i+ 1 - x i) ^r(«i+l + «*)/ 2 


=- K«i+1 - Ui) - p( x i+i - x i). 


+ 


(*£i+ 1 2'i)('5i+l T Sj)/2 


0 


n—l,n 

, (56) 


where 


T 
± r 


i-J rp 

a + u/L r 


(57) 


Observe that we have p = V{ui + \ — Ui)/(xi+ 1 — x) in the pseudo steady state. 
The cell- residual is then distributed to the nodes in all cells to yield the residual 
equation at node i: 

0 = + B- R t R ), (58) 

where the pseudo-time derivatives have been dropped. 
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We note that the Jacobian of this nonlinear advection-diffusion system 
is slightly different than that of the linear problem. The Jacobian blocks, 
however, are function of derivative of the distribution matrices B~ and B + : 


Ji—l 


Ji 


Jr 


i + 1 


1 d (B+<S> L ) 
hi dUi. _! ’ 

i ( a (Bt^) 

h ^ au, aut j ' 

i a (b- r * r ) 

hi dJJ i - |_i 


(59) 

(60) 
(61) 


In the special case a/u = constant, which is the only case considered in this pa- 
per, the distribution matrices reduce entirely to constant matrices. Therefore, 
we can entirely ignore the derivatives of B~ and Bf in deriving the Jaco- 
bian. As a result, the Jacobians are given formally as in Eq. (26). This makes 
the implementation of both linear and nonlinear systems extremely easy and 
straightforward. Effects of nonlinearity arise only in the derivatives of the left 
and right cell-residuals: 


d$ R 

~dU) 


Ori — Hr 


a 


U i+ 1 - Ui du u(u i+ 1 - Ui ) 


T 
± r 


dui 


+ 


2A t 
ph R 


-1 


k 


da 


+ 


1 du 


du.; L r dm 


v 


T 
± r 


h 


R 


2 T r J 
(62) 


d$ L 

~dU) 


~Ri 


hr 


a 


k 


u , 


'2 At 

-u ^ i du u(ui- Ui-i) —ph L ( da 


T 

± r 


du, 


+ 


- 


+ 


1 du 


V dui L r du. 


u 

+ = 
Tr 


hr 


2 T r J 
(63) 

where the average values indicated by the over bar are understood to be taken 
over the corresponding cell. Note that these derivatives depend on the solution 
at the current iteration, and therefore the superscript k has been introduced. 
The same applies to the distribution matrices also. Obviously, the nonlinear 
system approaches the linear system in the event that both the advection and 
the diffusion coefficients are constant. Also, the boundary conditions can be 
implemented as described in Section 3.3 for the linear case. Therefore, our 
hyperbolic nonlinear time dependent advection-diffusion system scheme can 
be implemented for linear, nonlinear, steady and unsteady problems without 
loss of accuracy. For simplify the discussion we only showed Jacobian for non- 
linear cases with a/u = constant but the scheme works for general nonlinear 
problems. 
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6 Results 


In this section we present the results in three categories: 1) Steady advection- 
diffusion equation for high Reynolds (or Peclct) number applications, 2) Un- 
steady linear advection-diffusion, and 3) Unsteady nonlinear advection-diffusion 
problems. We then present the order of accuracy results in the last subsec- 
tion. Note that time-accurate computations are started by BDF1 in the first 
time step with extremely small time step (e.g. t = 10” 8 ), and then by BDF2 
thereafter with much larger time step (e.g. t = 0.001 — 0.5) for all unsteady 
cases. This will ensure that the order of accuracy stays second order through 
all times. We remark that explicit time stepping is not available for time- 
accurate computations with the hyperbolic system method. 

6.1 Steady Linear Advection-Diffusion (High Re Appli- 
cations) 

Consider the advection-diffusion equation in x G (0, 1) with w(0) = 0 and 
■u(l) = 1 boundary conditions: 

d t u + ad x u = vd xx u + s(x), (64) 

where 

7 r 

s(x) = — [acos(7nr) + nu sin(7nc)], Re = a/v. (65) 

Re 

This is a boundary layer problem with a non-trivial steady state solution in 
the diffusion limit as a result of the source term addition [2] . This equation 
develops a very narrow boundary layer near the right boundary (x — 1) where 
the advection becomes dominant. The exact steady state solution to this 
problem is given by 

p—Re p(x—l)Re i 

u exact {x) = e _ Re _ + — sin(Tnr). (66) 

We chose Re values ranging from 1 to 10 6 and solved the equation on nonuni- 
form grid sizes up to 30000 nodes. We remark that the high-Re cases re- 
quired extremely fine grids to meet the well-known requirement on the mesh 
Reynolds-number [2], If desired, the computations can be performed on sub- 
stantially coarser grids with more aggressive and customized grid stretching 
as well as with some nonlinear mechanism to prevent numerical oscillations. 
However, we simply refined the grid to meet the mesh Reynolds-number re- 
quirement because our method is powerful enough to solve the problem very 
efficiently (i.e., less than 10 Newton iterations) even on such dense grids. The 
ability to efficiently solve the problem on highly refined grids is a great ad- 
vantage of these schemes. Shown in Fig. 5 are the comparisons between the 
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Figure 5: Residual-Distribution scheme on hyperbolic advection-diffusion sys- 
tem for boundary layer equation. 


computational and the exact solutions. In this comparison, the exact analytic 
solution is represented with + (plus) symbols while the numerical results are 
shown with solid lines (denoted Hyp-ADE). The exact solutions are plotted at 
nodes, and therefore they also show the computational grids. The solutions 
were obtained with the Newton-GS method as described in Section 3.2. For 
this problems, within each Newton iteration, the GS relaxation were conducted 
until three orders of magnitude reduction is achieved in the linear solver. The 
computations were continued until the residuals of both the solution and the 
solution gradient were reduced by eight orders of magnitude. Table 1 gives 
the convergence data obtained using the implicit formulation of the hyper- 
bolic advection-diffusion scheme. The linear dependency of the iterations on 

Table 1: Boundary layer problem (Convergence criteria: Residuals < ICR 8 .) 


logio Re 

Number of nodes 

GS relaxations/Newton iteration 

Newton iteration 

0 

10 

17 

7 

0 

100 

324 

7 

0 

200 

631 

7 

0 

300 

967 

7 

1 

300 

549 

7 

2 

300 

131 

7 

3 

300 

28 

7 

4 

300 

38 

6 

5 

3000 

60 

6 

6 

30000 

50 

5 


the grid size is demonstrated at Re = 1, which is a consequence of solving the 
advection-diffusion equation as a hyperbolic system for all Reynolds numbers 
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through the diffusion limit. This is remarkable because the linear convergence 
is retained for any irregular grid in any dimensions (N is approximately the 
number of nodes in each coordinate direction in two and three dimensions) 
as demonstrated in Refs. [1-5]. In conventional methods, e.g., such as the 
alternating direction implicit method with Thomas’ algorithm, the linear con- 
vergence can be achieved only on Cartesian grids (one-dimensional grid must 
exist in each coordinate direction). The O(N) convergence on general prob- 
lems and arbitrary grid systems makes the current scheme extremely attractive 
because it leads to orders of magnitude faster convergence in comparison with 
conventional methods whose convergence is typically 0(N 2 ). 

Table 1 shows for the grid of 300 nodes that the number of GS relaxations 
decreases significantly as Re increases. In Ref. [2] , a formula for L r was derived 
that minimizes the condition number of the hyperbolic advection-diffusion 
system and a nearly uniform convergence was obtained for an explicit pseudo- 
time marching scheme. We tested the formula in our method, but found that it 
did not change the significant decrease although it did reduce the number of GS 
relaxations approximately by 10 %. A different approach may be necessary to 
derive a formula to achieve Re-independent convergence of the GS relaxation 
in the implicit solver. Further study on Re-independent convergence is left 
as future work. Nevertheless, it should be noted that the number of Newton 
iterations is nearly Re-independent and very small: only 5 to 7 iterations to 
reduce the residual by eight orders of magnitude. Considering the fact that the 
cost of one GS relaxation is significantly cheaper than one Newton iteration, 
we find that the developed solver is still a very powerful steady solver. It thus 
serves as a very efficient solver for the unsteady residual equations as shown 
in the following sections. 

Finally, we remark that the high-Re cases required extremely fine grids to 
meet the well-known requirement on the mesh Reynolds-number as described 
in Ref. [2], If desired, the computations can be performed on substantially 
coarser grids with more aggressive grid stretching. However, we simply refined 
the grid to meet the mesh Reynolds-number requirement because our method 
is powerful enough to solve the problem very efficiently (i.e., 5 to 7 Newton 
iterations) even on such dense grids. The ability to efficiently solve the problem 
on highly refined grids is a great advantage. 


6.2 Unsteady Linear Advection-Diffusion 

We present unsteady computations for the two types of boundary conditions; 
Dirichlet (or Neumann) and periodic. We again note that for hyperbolic sys- 
tem scheme the implementation of Dirichlet and Neumann boundary condi- 
tions are similar. 
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Consider the time-dependent advection-diffusion equation in x G (0, 1) 


d t u + ad x u = vd xx u, (67) 

with the following time-dependent boundary condition (Dirichlet): 

u(0,t) = 0, (68) 

u(l,t) = U cos (cot), (69) 


where U is the amplitude of the oscillation and u is the frequency of the os- 
cillation on the right boundary. This problem has the following exact analytic 
solution: 


u 


exact 


(x, t ) = Real 


0 Xix 


0 \2X 


o X \ 


0X2 


-Ue 


iu)t 


a ± a / a 2 + Aiujis 

’ 1,2 - 2u ’ ^ ' 


where i = y/^1. We solved this time-dependent advection-diffusion equation 
with the following first-order hyperbolic system: 


d T u + ad x u = 

„ 3 Au n - u n ~ l 

uo r p — u -t , 

F 2A t 2 At 

( 71 ) 

d T p = 

(d x u — p)/T r . 

(72) 


The solution evolution was started with the exact initial solution given as 

( „\lX _ „\2X \ 

— in ) ’ (73) 

and was advanced in physical time after the time-dependent residual was re- 
duced by two orders of magnitude at each physical time step. We tested the 
scheme on a number of nonuniform grid systems ranging from N = 100... 10000 
nodes. Shown in Fig. 6 are the solution and the solution gradient obtained 
on the coarsest mesh (N=100) for U — 1, v — a — 1, and u = 77t/2. We 
show the results on the coarsest mesh just to demonstrate the accuracy of the 
scheme. Similar results were obtained on finer grid systems. The results are 
then compared with the analytic solution. Here the numerical results (denoted 
Hyp-ADE) are over plotted with the exact solutions, meaning excellent agree- 
ment. We will show the residual values in the accuracy verification section. 

We also obtained convergence data for different grid systems over 100 time 
steps. These are tabulated in Table 2, which shows a linear increase in number 
of iterations with increase of the grid system. This is, again, a consequence of 
solving the advection-diffusion equation as a hyperbolic system, where there 
are only first-derivatives and therefore the number of iterations is 0(N ), not 
0(N' 2 ). It is also shown in Table 2 that only four Newton iterations are 
required for all grids considered. 
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Figure 6: Time-dependent advection-diffusion with oscillatory boundary con- 
dition u(0,t) = 0, u(l,t) = cos(-yt) on N = 100 irregular nodes. Time step = 
0.01s. Left figures at t — 0.1s, right figures at t — 0.4s. 


Table 2: ffyperbolic time-dependent linear advection-diffusion problem with 
BDF2. Average data over 100 time steps are given. (Convergence criteria: 
Residuals < 10~ 2 ) 


Number of nodes 

GS relaxations/Newton iteration 

Newton iterations 

100 

127 

4 

300 

370 

4 

500 

607 

4 

1000 

1201 

4 
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For the periodic boundary condition demonstration, consider the time- 
dependent advection-diffusion equation in x G (0,1) given as in Eq. (67) with 
the following initial condition: 

u(x,t — 0) = sin(ra), (74) 

where k is an arbitrary constant. The exact solution to this problem with a 
periodic boundary condition is: 

U exact (x, t ) = e - * 2 " 4 sin(/c(x — at)). (75) 

We solved this problem with the same first-order hyperbolic advection-diffusions 
equation given as in Eqs. (71) and (72). For each physical time, we reduced 
the residuals by two orders of magnitude before advancing in time. During 
each time step, we also relaxed the linear system using GS relaxations until 
two orders of magnitude reduction in the linear system residuals was achieved. 
Shown in Fig. 7 are the solution and solution gradient on N — 25 uniform grid 





Figure 7: Time-dependent linear advection-diffusion problem(a = 1, v = 0.03) 
with periodic boundary condition on N = 25 uniform nodes. Time step = 
0.01s. From left to right, solutions at t — 0.1s, 0.5s, 1.0s. 

for a = 1, v = 0.03, and k = 2n. The results (Hyp-ADE) are over plotted 
with the exact solution, indicating excellent accuracy on such a coarse grid. 
Similar results were also obtained for larger grid systems. 

We examined the convergence rate of this problem on several grid systems. 
Given in Table 3, are the average numbers of GS relaxations per Newton itera- 
tion obtained over 100 time steps. Clearly the convergence rate is of O(N), not 
0(N 2 ) as typical for numerical methods for the advection-diffusion equation. 
Observe that for most grid systems only one or two Newton iterations were 
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Table 3: Hyperbolic time-dependent linear advection-diffusion problem with 
periodic BC (a = l,u = 0.03) Average data over 100 time steps are 
given. (Convergence criteria: Residuals < 10~ 2 ) 


Number of nodes 

GS relaxations/Newton iteration 

Newton iterations 

25 

4 

3 

100 

11 

3 

300 

33 

2 

500 

55 

2 

1000 

116 

2 


sufficient to obtain accurate solutions. We also note that we used At = 0.01- 
sec for all the grid systems; the BDF schemes are unconditionally stable. The 
time step is orders-of-magnitude larger than that required for conventional 
explicit schemes, which is limited by 0(h 2 ). Of course, conventional implicit 
schemes also allow unconditionally large time steps, but it requires 0(N 2 ) 
convergence in an iterative linear solver and potentially a much larger num- 
ber of outer iterations as well if the exact linearization is not possible and 
Newton’s method cannot be constructed. Hence, the method developed here 
has two major advantages over conventional methods: the exact linearization 
(Newton’s method) and O(N) iterative convergence in the linear solver. The 
latter advantage can be potentially huge with increase of the grid system as 
the speed-up factor is O(N) and thus grows for finer grids. For example, for 
N = 1000, our scheme is at least four orders of magnitude faster than the 
conventional schemes, which is remarkable. 


6.3 Unsteady Nonlinear Advection-Diffusion 

Consider the time-dependent nonlinear viscous Burgers equation with a time- 
dependent source term: 

d t u + d x f = d x {vu x ) + g(x,t), x e (0, 1), (76) 

where / = u 2 / 2, v — u, and 

g(x,t) = u e t + ^((u e ) 2 ) x - «) 2 -u e u e xx . (77) 

The source term has been generated by the following function: 


u e (x, t ) = Real 


SmHx ^) Ue ^) +CiC>1 , 


smn \/i(jj v) 
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so that it is the exact solution to Eq. (76). The boundary conditions are 


u(0,t) = C, (79) 

u(l,t) = C + Ucos(cjt), (80) 

where uu is the frequency of the oscillation on the right boundary, and U is the 
amplitude of the oscillation. The constant C must be greater than 1 to keep 
the diffusion coefficient positive. 

We solved this time-dependent nonlinear advection-diffusion equation with 
the following equivalent first-order hyperbolic system: 

3 _ u n-i 

d T u + d x ( u 2 ) = d x p - —u + — + g(x, t), (81) 

^ d T p = (d x u-p/u). (82) 

Shown in Fig. 8 are the results obtained on N = 100 nonuniform nodes for 
U = 1, uj = 77t/2, and C — 2. The scheme produces very accurate time 



Figure 8: Time-dependent nonlinear advection-diffusion (a = v = u) with 
oscillatory boundary condition w(0,f) = 2, u(l,t) = 2 + cos(^-t) on N = 100 
nonuniform nodes. Time step = 0.01s. From left to right, solutions at t = 
0.1s, 2.5s, 10.0s. 

evolutions of the solution (u) and the gradient (p/v) on nonuniform grids (see 
also Fig. 9.) 

The 0(N ) convergence rate was once again achieved for the time-dependent 
nonlinear hyperbolic advection-diffusion system. This is given in Table 4, 
where the average number of iterations were obtained over 1000 time steps 
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Figure 9: Time-dependent nonlinear advection-diffusion (a = v = u) with 
oscillatory boundary condition u(0,t) = 2, u(l,t) = 2 + cos(A-t) on N = 100 
nonuniform nodes. Time step = 0.01s 


(over 17 periods). Note also that the method converged at four Newton iter- 
ations for all grids considered. Again, a constant time step (At = 0.01-sec) 
was used for all the grid systems. And the method is substantially more ef- 
ficient here also than conventional explicit /implicit schemes as discussed for 
the linear case in Section 6.2. 

Table 4: Hyperbolic time-dependent nonlinear advection-diffusion problem 
with BDF2. Average data over 100 time steps are given. (Convergence criteria: 
Residuals < 10” 2 ) 


Number of nodes 

GS rclaxations/Newton iteration 

Newton iterations 

100 

88 

4 

300 

263 

4 

500 

440 

4 

1000 

882 

4 


6.4 Accuracy Verification 

To demonstrate the overall order of accuracy of our time-dependent hyperbolic 
advection-diffusion scheme, we computed the L x = YliLi (U^ xact — Ui)/N as well 
as Lqo = Max(Uf xact — Ui) norms of errors on a series of nonuniform grids of 
N nodes, and time step increments At. For each case, we terminate the solver 
when the residuals were dropped by two orders of magnitude. We remark 
that in some practical applications two orders of magnitude reduction in the 
residuals may not be sufficient to achieve the appropriate order of accuracy. 
But for the cases considered here, terminating at even nine orders of magnitude 
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resulted in similar order of accuracy. The spatial accuracy was obtained by 
comparing several nonuniform grid systems on a fixed At, while the temporal 
accuracy was obtained by comparing various At on a fixed nonuniform grid 
system. We also included boundary nodes in both Li and norms. For 
unsteady cases, the L\ and L ^ norms were evaluated at t — 1.0-sec. 


2 -3 


Figure 10: error obtained for time-dependent linear hyperbolic advection- 

diffusion system scheme. Left: Spatial accuracy, Right: Temporal accuracy. 

Figures 10 and 11 show the L ^ error convergence results for both linear 
and nonlinear time-dependent hyperbolic advection diffusion system, where h 
is the representative mesh spacing defined by h = / (N — 1). For discussion 
purposes, we only present the accuracy plots for Re = 1 results; similar results 
were obtained for other Reynolds numbers. These results show that our scheme 
is second-order accurate for all the variables and the gradients (note that again 
that our flux integral is exact) at all the grid nodes including the boundary 
nodes (see also Tables 5, 6, 7, and 8). We note that it is natural for the error 
convergence in the temporal space to deteriorate when the 0(At 2 ) « 0(h 2 ) 
as well as when At is too large for the time discretization to be in asymptotic 
range; i.e. a lower order of accuracy for coarse discretization in time is typical. 
These limits are also included in the temporal accuracy plot for completeness. 

Table 5: Spatial accuracy for the linear advection-diffusion problem with oscil- 
latory boundary condition on several nonuniform grid systems (At = 0.001s). 




Number of nodes 

L\ error of u 

Loo error of u 

Order 

Li error of p 

Loo error of p 

Order 

25 

9.65E-04 

3.44E-03 


6.58E-03 

1.71E-02 


50 

2.28E-05 

8.98E-04 

1.94 

1.08E-03 

4.22E-03 

2.02 

100 

6.10E-05 

2.13E-04 

2.08 

3.08E-04 

9.95E-04 

2.08 

200 

1.59E-06 

5.04E-05 

2.08 

8.76E-05 

2.35E-04 

2.08 

300 

7.58E-06 

2.06E-05 

2.21 

4.87E-05 

9.53E-05 

2.23 
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Figure 11: error obtained for time-dependent nonlinear hyperbolic 

advection-diffusion system scheme. Left: Spatial accuracy, Right: Temporal 
accuracy. 


Table 6: Temporal accuracy for the linear advection-diffusion problem with 
oscillatory boundary condition on a nonuniform grid (N=100) with huax. = 


0.045. 


At 

Li error of u 

Loo error of u 

Order 

Li error of p 

Loo error of p 

Order 

0.500 

3.01E-01 

6.63E-01 


2.46E+00 

4.09E+00 


0.100 

2.32E-02 

7.55E-02 

1.35 

1.25E-01 

2.32E-01 

1.78 

0.050 

4.37E-03 

1.72E-02 

2.13 

3.57E-02 

6.28E-02 

1.88 

0.025 

9.66E-04 

3.16E-03 

2.44 

1.16E-02 

2.26E-02 

1.47 

0.010 

1.48E-04 

2.99E-04 

2.57 

2.21E-03 

4.63E-03 

1.73 

0.001 

6.02E-05 

2.13E-04 

0.15 

3.05E-04 

9.95E-04 

0.67 


Table 7: Spatial accuracy for the nonlinear advection-diffusion problem with 
oscillatory boundary condition on several nonuniform grid systems (At = 
0.001s). 


Number of nodes 

Li error of u 

Loo error of u 

Order 

Li error of p 

Lqo error of p 

Order 

10 

25 

1.83E-03 

4.36E-04 

8.41E-03 

1.50E-03 

1.88 

1.20E-02 

2.79E-03 

4.46E-02 

7.29E-03 

1.98 

50 

9.00E-05 

3.74E-04 

2.00 

5.48E-04 

1.78E-03 

2.03 

100 

2.80E-05 

9.39E-05 

1.99 

1.41E-04 

4.36E-04 

2.03 

200 

8.19E-06 

2.33E-05 

2.01 

4.81E-05 

1.07E-04 

2.03 
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Table 8: Temporal accuracy for the nonlinear advection-diffusion problem 
with oscillatory boundary condition on a nonuniform grid system (N=100) 
with JiMax. = 0.045. 


At 

Li error of u 

Loo error of u 

Order 

Li error of p 

Loo error of p 

Order 

0.200 

6.46E-02 

1.55E-01 


4.97E-01 

6.83E-01 


0.100 

1.52E-02 

4.67E-02 

1.73 

9.00E-02 

1.54E-01 

2.15 

0.050 

2.03E-03 

7.40E-03 

2.66 

1.99E-02 

3.45E-02 

2.16 

0.020 

3.52E-04 

7.32E-04 

2.52 

5.01E-03 

1.03E-02 

1.32 

0.010 

1.34E-04 

3.20E-04 

1.19 

1.70E-03 

3.06E-03 

1.75 

0.005 

6.68E-05 

1.41E-04 

1.18 

6.20E-04 

9.90E-04 

1.63 

0.001 

2.80E-05 

9.39E-05 

0.25 

1.41E-04 

4.36E-04 

0.51 


7 Conclusions 

In this paper, we have extended the first-order hyperbolic system method to 
time-accurate computations. The time integration scheme is implicit with the 
second-order backward-difference formula, and the resulting system of time- 
dependent residual equations is solved by a steady solver developed based on 
the hyperbolic method. The steady solver is Newton’s method constructed 
from a compact upwind residual-distribution scheme for linear and nonlinear 
advection-diffusion equations. The construction of the numerical scheme was 
greatly simplified for the nonlinear equation by the preconditioned conservative 
formulation. We have demonstrated that the developed time-accurate schemes 
are capable, at every physical time step, of producing second-order accurate 
solution and gradient on highly refined nonuniform grids by less than five 
Newton iterations with the linear convergence of the block GS relaxations for 
all test cases considered. Thus, the notable features of the hyperbolic method 
have been shown to carry over to time-dependent problems. 

The developed methodology based on the second-order backward Euler 
time integration can be extended to practical one-dimensional problems [7-11], 
bringing significant improvements in both accuracy and efficiency. Extension 
to higher-order accuracy is also possible by improving the residual evaluation 
and higher-order backward-difference formulas. It is relatively simple in one 
dimension because the integration of the flux balance term is exact in one 
dimension as mentioned in Section 3.1. In higher dimensions, although the 
construction of time-accurate schemes is as simple as presented in the paper, 
it is known that a high-order quadrature would be required for the flux balance 
integration in order to achieve the equal order of spatial accuracy in the solu- 
tion and the gradient [2], Although not considered in the present paper, the 
construction of non-oscillatory schemes is an important area of development 
for problems with discontinuities such as shock waves. 
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